import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('dark_background')
Use the starter program ising.py (or its Matlab version ising.m / calculate_spin.m) or your own equivalent routines to simulate the Ising ferromagnet in 2D. Do production runs at temperatures $T$ sampled in regular increments both above and below $T_C\approx2.27$ (measured in units of J/kB), for a few different values of the magnetic field $H$ (measured in units of J/$\mu$). Make sure to calculate for $H=0$, as well as positive and negative $H$. For initial conditions, you can try, e.g., having all spins up. Observe and note how the magnetization per spin ($m$) and energy per spin ($E$) converge towards their equilibrium values, as well as the qualitative features of the equilibrium configurations that may be reached. For example, check how the rate of convergence depends on the temperature and/or the field.
Using code from my GitHub repo on the subject (Stripped down version found below):
# After running `.\IsingModel`
a = np.load('results.npz')['isingmodel']
df = pd.DataFrame(a, columns=['Temperature', 'MagneticField', 'Energy', 'Magnetization', 'SpecificHeatCapacity', 'MagneticSusceptibility'])
df
| Temperature | MagneticField | Energy | Magnetization | SpecificHeatCapacity | MagneticSusceptibility | |
|---|---|---|---|---|---|---|
| 0 | 1.500000 | 0.0 | -3.899840 | 0.985893 | 0.783993 | 0.027264 |
| 1 | 1.502002 | 0.0 | -3.904960 | 0.986933 | 0.801392 | 0.026147 |
| 2 | 1.504004 | 0.0 | -3.899058 | -0.986151 | 0.798130 | 0.026580 |
| 3 | 1.506006 | 0.0 | -3.902969 | -0.986462 | 0.750168 | 0.027782 |
| 4 | 1.508008 | 0.0 | -3.902329 | 0.986542 | 0.721839 | 0.024825 |
| ... | ... | ... | ... | ... | ... | ... |
| 995 | 3.491992 | 0.0 | -1.344462 | 0.013209 | 1.162247 | 2.146309 |
| 996 | 3.493994 | 0.0 | -1.321209 | 0.023609 | 0.987977 | 1.791802 |
| 997 | 3.495996 | 0.0 | -1.305138 | -0.010960 | 1.059259 | 1.819019 |
| 998 | 3.497998 | 0.0 | -1.305351 | -0.018427 | 1.052455 | 1.599211 |
| 999 | 3.500000 | 0.0 | -1.298987 | -0.000542 | 0.923059 | 1.483989 |
1000 rows × 6 columns
fig, ax = plt.subplots(2,2, figsize=(16, 16), dpi=300, constrained_layout=True)
fig.suptitle('Ising Model (Metropolis-Hastings Algorithm)')
df.plot.scatter('Temperature', 'Energy', ax=ax[0,0], s=0.1)
df.plot.scatter('Temperature', 'Magnetization', ax=ax[0,1], s=0.1)
df.plot.scatter('Temperature', 'SpecificHeatCapacity', ax=ax[1,0], s=0.1)
df.plot.scatter('Temperature', 'MagneticSusceptibility', ax=ax[1,1], s=0.1)
plt.show()
# After running `.\IsingModel --nH 100 --nT 100 --Hi -10 --Hf 10`
a = np.load('results_3d.npz')['isingmodel']
df = pd.DataFrame(a, columns=['Temperature', 'MagneticField', 'Energy', 'Magnetization', 'SpecificHeatCapacity', 'MagneticSusceptibility'])
df
| Temperature | MagneticField | Energy | Magnetization | SpecificHeatCapacity | MagneticSusceptibility | |
|---|---|---|---|---|---|---|
| 0 | 1.5 | -10.000000 | -14.000000 | 1.000000 | 0.000000e+00 | 0.000000 |
| 1 | 1.5 | -9.797980 | -13.797980 | 1.000000 | 2.359405e-24 | 0.000000 |
| 2 | 1.5 | -9.595960 | -13.595960 | 1.000000 | 5.308661e-24 | 0.000000 |
| 3 | 1.5 | -9.393939 | -13.393939 | 1.000000 | 6.903384e-26 | 0.000000 |
| 4 | 1.5 | -9.191919 | -13.191919 | 1.000000 | 2.680065e-24 | 0.000000 |
| ... | ... | ... | ... | ... | ... | ... |
| 9995 | 3.5 | 9.191919 | -13.173123 | -0.998907 | 4.884310e-02 | 0.000578 |
| 9996 | 3.5 | 9.393939 | -13.378942 | -0.999138 | 4.197148e-02 | 0.000486 |
| 9997 | 3.5 | 9.595960 | -13.582196 | -0.999218 | 3.875758e-02 | 0.000438 |
| 9998 | 3.5 | 9.797980 | -13.779470 | -0.998960 | 5.576773e-02 | 0.000616 |
| 9999 | 3.5 | 10.000000 | -13.989120 | -0.999396 | 3.168047e-02 | 0.000342 |
10000 rows × 6 columns
fig = plt.figure(figsize=(16,16), dpi=300, constrained_layout=True)
fig.suptitle('Ising Model (Metropolis-Hastings Algorithm)')
ax1 = fig.add_subplot(221, projection='3d')
ax1.view_init(30, 120)
ax1.plot_trisurf(df.MagneticField, df.Temperature, df.Energy, cmap=plt.get_cmap('jet'))
ax1.set_xlabel('Magnetic Field')
ax1.set_ylabel('Temperature')
ax1.set_zlabel('Energy')
ax2 = fig.add_subplot(222, projection='3d')
ax2.view_init(30, 70)
ax2.plot_trisurf(df.MagneticField, df.Temperature, df.Magnetization, cmap=plt.get_cmap('jet'))
ax2.set_xlabel('Magnetic Field')
ax2.set_ylabel('Temperature')
ax2.set_zlabel('Magnetization')
ax3 = fig.add_subplot(223, projection='3d')
ax3.view_init(30, -110)
ax3.plot_trisurf(df.MagneticField, df.Temperature, df.SpecificHeatCapacity, cmap=plt.get_cmap('jet'))
ax3.set_xlabel('Magnetic Field')
ax3.set_ylabel('Temperature')
ax3.set_zlabel('Specific Heat Capacity')
ax4 = fig.add_subplot(224, projection='3d')
ax4.view_init(30, -110)
ax4.plot_trisurf(df.MagneticField, df.Temperature, df.MagneticSusceptibility, cmap=plt.get_cmap('jet'))
ax4.set_xlabel('Magnetic Field')
ax4.set_ylabel('Temperature')
ax4.set_zlabel('Magnetic Susceptibility')
plt.show()